knitr::opts_chunk$set(echo = TRUE, warning = FALSE)
options(width = 100) # ширина текстового вывода
options(digits = 6) # число знаков после запятой в выводе
library(lubridate)
##
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
##
## date
library(ggrepel)
## Loading required package: ggplot2
library(readr)
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.2.1 --
## v tibble 2.1.3 v dplyr 0.8.3
## v tidyr 1.0.0 v stringr 1.4.0
## v purrr 0.3.2 v forcats 0.4.0
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x lubridate::as.difftime() masks base::as.difftime()
## x lubridate::date() masks base::date()
## x dplyr::filter() masks stats::filter()
## x lubridate::intersect() masks base::intersect()
## x dplyr::lag() masks stats::lag()
## x lubridate::setdiff() masks base::setdiff()
## x lubridate::union() masks base::union()
library(ggplot2)
library(ggfortify)
library(modelr)
library(broom)
##
## Attaching package: 'broom'
## The following object is masked from 'package:modelr':
##
## bootstrap
library(GGally)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
##
## Attaching package: 'GGally'
## The following object is masked from 'package:dplyr':
##
## nasa
library(forcats)
library(tidyr)
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
library(dplyr)
library(memisc)
## Loading required package: lattice
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
##
## Attaching package: 'memisc'
## The following object is masked from 'package:car':
##
## recode
## The following objects are masked from 'package:dplyr':
##
## collect, recode, rename, syms
## The following object is masked from 'package:purrr':
##
## %@%
## The following object is masked from 'package:ggplot2':
##
## syms
## The following object is masked from 'package:lubridate':
##
## is.interval
## The following objects are masked from 'package:stats':
##
## contr.sum, contr.treatment, contrasts
## The following object is masked from 'package:base':
##
## as.array
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
datetime: дата и время аренды велосипеда;season: время года (1 - весна, 2 - лето, 3 - осень, 4- зима);holiday: выходной (1 - да, 0 - нет);workingday: рабочий день (1 - да, 0 -нет);weather: качественная характеристика погодных условий:temp: температура в цельсиях;atemp: температура в цельсиях (по ощущениям);humidity: относительная влажность;windspeed: скорость ветра;casual: число актов аренды велосипедов незарегестрированными пользователями;registered: число актов аренды велосипедов зарегестрированными пользователями;count - целевая переменная: общее число актов аренды;bike <- read_csv ("C:/Users/Юра/Desktop/KP-NIS-R/ПРОЕкт/Проек/bike.csv")
## Parsed with column specification:
## cols(
## datetime = col_datetime(format = ""),
## season = col_double(),
## holiday = col_double(),
## workingday = col_double(),
## weather = col_double(),
## temp = col_double(),
## atemp = col_double(),
## humidity = col_double(),
## windspeed = col_double(),
## casual = col_double(),
## registered = col_double(),
## count = col_double()
## )
Вывод: была загружена база данных по арендам велосипедов. Она содержит в себе 10886 наблюдений, которые представляют собой акты аренды велосипедов с предшествующими ими внешними условями, такими как погодные условия, сезон, календарный день и т.д - над ней будет проводиться анализ и поиск подходящей модели прогнозирования. Нам необходимо подготовить данные к анализу и разбить базу данных на обучающую и тестовую выборки случайным образом, и именно к тестовой выборке будет применяться выявленная на обучающей выборке лучшая модель для прогнозирования значений актов аренды велосипедов.
Удалим переменные “registered” и “casual” в нашей выборке, потому что они являются лишь слагаемыми нашей целей переменной “count”, поэтому не несут в себе большой значимости и могут изменить модели корреляции по причиние их сильной зависимости от целевой переменной:
bike <- bike[,-(10:11)]
Удалим времена года и выделим из дат месяца в обучающей и тестовой выборках, а также время, чтобы повысить точность нашего анализа, который может быть основан на месяцах, а не просто временах года:
## Удаление времён года
bike <- bike[,-2]
## Создание отдельных столбцов с временем и датой в выборках
bike$month <- as.Date(bike$datetime)
bike$time <- format(as.POSIXct(bike$datetime) ,format = "%H:%M:%S")
## Удаление объединённого столбца в выборках
bike <- bike[,-1]
## Выделение месяцев из дат
bike$month <- format(bike$month, "%m")
## Выделение часов из времени
bike$time <- substr(bike$time, start = 1, stop = 2)
bike$month <- as.numeric(as.character(bike$month))
bike$time <- as.numeric(as.character(bike$time))
bike <- mutate (bike,
daytime = ifelse(time >= 0, ifelse(time<=6, 1,ifelse(time<=12, 2,ifelse (time <=17,3, 4))), 0))
## Удаление ненужного столбца с часами
bike <- bike[,-10]
holiday: выходной (1 - да, 0 - нет);workingday: рабочий день (1 - да, 0 -нет);weather: качественная характеристика погодных условий:temp: температура в цельсиях;atemp: температура в цельсиях (по ощущениям);humidity: относительная влажность;windspeed: скорость ветра;month: месяц года (от 1 до 12);daytime: время суток: 1: ночь; 2: утро; 3: день; 4: вечер;count - целевая переменная: общее число актов аренды;sample <- floor(0.8*nrow(bike))
set.seed(10886)
sample_ind <- sample(seq_len(nrow(bike)),size = sample)
bike_train <- bike[sample_ind,]
bike_test <- bike[-sample_ind,]
Прежде чем начать моделирование, необходимо выяснить: как сильно на корреляцию влияют удаленные столбцы registered и casual.Так как целевая переменная count = registered+casual, очевидно, что значение двух нецелевых переменных влияют напрямую на целевую. Задание по прогнозированию: прогноз переменной count. Поэтому, при анализе и моделировании необходимо исключить registered и casual. Однако, чтобы не потерять скрытые зависимости, которые несут в себе удаленные столбцы, проведем анализ этих столбцов на зависимости. В лучшем исходе registered и casual должны показывать одинаковые зависимости и тренды, что позволи потсроить модель без них и без ухудшения самой модели.
Загрузим изначальную выборку с переменными “registered” и “casual”. Сделаем отношение в 0.75, которое означает, что если в одном наблюдении более 75% арендаторов велосипедов были зарегистрированными, то считается, что во время этого наблюдения были только зарегистрированные кленты, если менее - то только незарегистрированные. Такая большая доля выбрана потому, что число незарегистрированных пользователей почти всегда незначительна.
bike_orig <- read_csv ("C:/Users/Юра/Desktop/KP-NIS-R/ПРОЕкт/Проек/bike.csv")
## Parsed with column specification:
## cols(
## datetime = col_datetime(format = ""),
## season = col_double(),
## holiday = col_double(),
## workingday = col_double(),
## weather = col_double(),
## temp = col_double(),
## atemp = col_double(),
## humidity = col_double(),
## windspeed = col_double(),
## casual = col_double(),
## registered = col_double(),
## count = col_double()
## )
bike_orig <- mutate(bike_orig, relation = ifelse(registered/count >= 0.75, "registered", "casual"))
graph1 <- bike_orig %>%
ggplot(aes(count, temp, color = relation)) +
geom_point(alpha = 0.1) +
geom_smooth(method = 'lm', se = F)+
labs(title = "1", colour = "Тип клиента", x = "Число актов аренды", y = "t воздуха")
graph2 <- bike_orig %>%
ggplot(aes(count, atemp, color = relation)) +
geom_point(alpha = 0.1) +
geom_smooth(method = 'lm', se = F)+
labs(title = "2", colour = "Тип клиента", x = "Число актов аренды", y = "t воздуха по ощущениям")
graph3 <- bike_orig %>%
ggplot(aes(count, weather, color = relation)) +
geom_point(alpha = 0.1) +
geom_smooth(method = 'lm', se = F)+
labs(title = "3", colour = "Тип клиента", x = "Число актов аренды", y = "Cостояние погоды")
graph4 <- bike_orig %>%
ggplot(aes(count, season, color = relation)) +
geom_point(alpha = 0.1) +
geom_smooth(method = 'lm', se = F)+
labs(title = "4", colour = "Тип клиента", x = "Число актов аренды", y = "Сезон")
grid.arrange(graph1, graph2, graph3, graph4, ncol=2)
Вывод:
1 график: можно заметить, что при высокой температуре воздуха существует тенденция для незарегистрированных пользователей для обращения к аренде велосипедов. Для зарегистрированных пользователей температура не так важна - их число актов аренды распределено равномерно;
2 график: закономерности практически идентичные графику №1;
3 график: можно заявить, что во время самой плохой степени погоды берётся крайне мало велосипедов, но вся их аренда производится только зарегистрированными пользователями. В остальных случаях, есть как зарегистрированные, так незарегистрированные пользователи (зарегистрированных по-прежнему больше);
4 график: можно заметить, что во время зимнего сезона арендуют велосипеды практически только зарегистрированные пользователи. Больше всего незарегистрированных пользователей арендуют велосипеды весной и летом;
graph5 <- bike_orig %>%
ggplot(aes(count, workingday, color = relation)) +
geom_point(alpha = 0.1) +
geom_smooth(method = 'lm', se = F)+
labs (title = "5", colour = "Тип клиента", x = "Число актов аренды", y = "Будние дни")
graph6 <- bike_orig %>%
ggplot(aes(count, holiday, color = relation)) +
geom_point(alpha = 0.1) +
geom_smooth(method = 'lm', se = F)+
labs(title = "6", colour = "Тип клиента", x = "Число актов аренды", y = "Праздничные дни")
graph7 <- bike_orig %>%
ggplot(aes(count, windspeed, color = relation)) +
geom_point(alpha = 0.1) +
geom_smooth(method = 'lm', se = F)+
labs(title = "7", colour = "Тип клиента", x = "Число актов аренды", y = "Скорость ветра")
graph8 <- bike_orig %>%
ggplot(aes(count, humidity, color = relation)) +
geom_point(alpha = 0.1) +
geom_smooth(method = 'lm', se = F)+
labs(title = "8", colour = "Тип клиента", x = "Число актов аренды", y = "Влажность")
grid.arrange(graph5, graph6, graph7, graph8, ncol=2)
Вывод:
5 график: можно заметить, что в выходные дни есть тенденция для незарегистрированных пользователей прибегать к аренде велосипедов, в то время как в буние дни берут велосипеды практически только зарегистрированные пользователи;
6 график: можно заметить, что в праздничные есть тенденция для незарегистрированных пользователей прибегать к аренде велосипедов, в то время как в обычные дни берут велосипеды практически только зарегистрированные пользователи;
7 график: нет практически никакой зависиомсти между зарегистрированностью пользователя и скоростью ветра. Все типы клиентов распределны случайным образом;
8 график: можно заметить, что при низкой влажности больше незарегистрированных пользователей берёт велосипеды в аренду. Тем не менее, эта зависимость выражена слабо - можно сказать, что влажность тоже не влияет на тенднецию типр=ов клиентов к аренде велосипедов;
Вывод: анализ зависимостей показал значительные несовпадения зависимостей переменных с рабочими днями, сезоном (месяцы) и влажность. Используя эти переменные при моделировании стоит быть аккуратными, т.к. их применение может исказить реальную ситуацию, что даст большую ошибку при прогнозирования.
summary (bike_train)
## holiday workingday weather temp atemp humidity
## Min. :0.0000 Min. :0.000 Min. :1.00 Min. : 0.82 Min. : 0.76 Min. : 0.0
## 1st Qu.:0.0000 1st Qu.:0.000 1st Qu.:1.00 1st Qu.:13.94 1st Qu.:16.66 1st Qu.: 47.0
## Median :0.0000 Median :1.000 Median :1.00 Median :20.50 Median :24.24 Median : 62.0
## Mean :0.0286 Mean :0.679 Mean :1.42 Mean :20.26 Mean :23.68 Mean : 61.9
## 3rd Qu.:0.0000 3rd Qu.:1.000 3rd Qu.:2.00 3rd Qu.:26.24 3rd Qu.:31.06 3rd Qu.: 77.0
## Max. :1.0000 Max. :1.000 Max. :4.00 Max. :41.00 Max. :45.45 Max. :100.0
## windspeed count month daytime
## Min. : 0.0 Min. : 1 Min. : 1.00 Min. :1.00
## 1st Qu.: 7.0 1st Qu.: 41 1st Qu.: 4.00 1st Qu.:1.00
## Median :13.0 Median :144 Median : 7.00 Median :2.00
## Mean :12.8 Mean :191 Mean : 6.53 Mean :2.42
## 3rd Qu.:17.0 3rd Qu.:283 3rd Qu.:10.00 3rd Qu.:4.00
## Max. :57.0 Max. :977 Max. :12.00 Max. :4.00
summary (bike_test)
## holiday workingday weather temp atemp humidity
## Min. :0.0000 Min. :0.00 Min. :1.0 Min. : 0.82 Min. : 1.51 Min. : 0
## 1st Qu.:0.0000 1st Qu.:0.00 1st Qu.:1.0 1st Qu.:13.94 1st Qu.:16.66 1st Qu.: 47
## Median :0.0000 Median :1.00 Median :1.0 Median :20.50 Median :24.24 Median : 62
## Mean :0.0285 Mean :0.69 Mean :1.4 Mean :20.13 Mean :23.54 Mean : 62
## 3rd Qu.:0.0000 3rd Qu.:1.00 3rd Qu.:2.0 3rd Qu.:26.24 3rd Qu.:31.06 3rd Qu.: 78
## Max. :1.0000 Max. :1.00 Max. :3.0 Max. :39.36 Max. :44.70 Max. :100
## windspeed count month daytime
## Min. : 0.0 Min. : 1 Min. : 1.00 Min. :1.00
## 1st Qu.: 7.0 1st Qu.: 46 1st Qu.: 3.00 1st Qu.:1.00
## Median :13.0 Median :148 Median : 7.00 Median :2.00
## Mean :12.8 Mean :193 Mean : 6.48 Mean :2.43
## 3rd Qu.:17.0 3rd Qu.:288 3rd Qu.: 9.00 3rd Qu.:4.00
## Max. :57.0 Max. :897 Max. :12.00 Max. :4.00
Вывод: статистический анализ переменных показал, что пустые ячейки отсутствуют, а также существуют возможные выбросы для переменных “windspeed” и “count” в связи большими отклонениями максимальных значений от 3 квартилей. Тем не менее, эти выбросы могут быть объяснены моделью и не являться для неё преградой, поэтому будет проведена отдельная работа с выбросами, которая выявит их значимость для работы моделей.
Также можно заметить, что в датасетах присутствуют 5 качественых переменных (holiday, workingday, weather, month и daytime) и 5 количественных переменных (temp, atemp, humidity, windspeed, count). Качественные переменные выражены количественными категориями, поэтому проведения работ по преобразованию их в “долгий” текстовый тип не нужно - они готовы к работе без дополнительных преобразований.Среди качественных переменных выбросов не было обнаружены, все их количественные предикторы лежат в опредеённых диапозонах, что исключает возможность опечатки при их заполнении.
Вывод: На диаграммах рассеяния видна зависимость между всеми параметрами. Коэффициенты корреляции ненулевые и статистически значимы. Наибольшее влияние на целевую переменную оказывают время суток,температура воздуха (по ощущениям тоже) и влажность. Наибольшая корреляция между переменными - temp и atemp, но это очевидно, ведь температура воздуха по ощущениям в наибольшей степени зависит от фактической. Влажность воздуха имеет достаточно большой коэффициент корреляции с типом погоды, что также может быть легко объяснено, в то время как скорость ветра в некоторой степени зависит от влажности, а значит и от погоды. Месяц года имеет хороший коэффициент корреляции с температурой воздуха и его влажностью, что также очевидно, ведь в разные времена года разные темепературные режимы и состояние воздуха в связи с их цикличностью. В остальных случаях, поверхностного анализа недостаточного - необходим более детальный анализ, либо проверка работы моделей с взаимосвязями переменных, что и будет проделано в рамках этого проекта.
model_1 <- lm(count ~ daytime, data = bike_train)
summary(model_1)
##
## Call:
## lm(formula = count ~ daytime, data = bike_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -293.8 -91.0 -51.7 68.3 739.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.38 4.10 6.68 2.5e-11 ***
## daytime 67.60 1.53 44.23 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 164 on 8706 degrees of freedom
## Multiple R-squared: 0.183, Adjusted R-squared: 0.183
## F-statistic: 1.96e+03 on 1 and 8706 DF, p-value: <2e-16
Вывод: время суток и число актов аренды показывали наибольшую корреляцию, поэтому было решено построить модель с такими зависимостями первой. Данная модель регрессии статистически значимая, так как p-value = 2e-16 (<0.05).
Проверка параметров модели:
glance(model_1)
## # A tibble: 1 x 11
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
## <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int>
## 1 0.183 0.183 164. 1956. 0 2 -56774. 113554. 1.14e5 2.35e8 8706
Вывод: модель смогла объяснить только 18.3% наблюдений, и её стандартная ошибка равна 164, что говорит о том, модель очень плохо справляется своими задачами и нуждается в замене альтернативами.
t1<- car::powerTransform(count ~ temp + atemp + daytime + humidity, data = bike_train)
coef(t1)
## Y1
## 0.273731
model_2 <- lm(count^0.274 ~ temp + atemp +I(atemp^2) + daytime +I(daytime^2) + humidity, data = bike_train)
summary(model_2)
##
## Call:
## lm(formula = count^0.274 ~ temp + atemp + I(atemp^2) + daytime +
## I(daytime^2) + humidity, data = bike_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1104 -0.5598 -0.0489 0.5207 2.5775
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.422615 0.079566 -17.88 <2e-16 ***
## temp 0.017194 0.006738 2.55 0.011 *
## atemp 0.085848 0.007242 11.85 <2e-16 ***
## I(atemp^2) -0.001291 0.000114 -11.30 <2e-16 ***
## daytime 3.304186 0.044506 74.24 <2e-16 ***
## I(daytime^2) -0.553026 0.008679 -63.72 <2e-16 ***
## humidity -0.006951 0.000481 -14.44 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.774 on 8701 degrees of freedom
## Multiple R-squared: 0.634, Adjusted R-squared: 0.634
## F-statistic: 2.52e+03 on 6 and 8701 DF, p-value: <2e-16
Вывод: Следующая модель - это модель с несколькими переменными, которые имеют наибольшие коэффициенты корреляции с целевой. Все, кроме одного коэффициента в модели, значимые, и сама она статистически значимая p-value < 2e-16 (<0.05).
Построение графиков остаточной диагностики:
autoplot (model_2)
Вывод:
Рассеяние точек говорит о наличии нелинейной связи между переменными (линейная связь накладывается на пунктирную прямую в нуле), но тем не менее, коэффициенты нелинейной модели, основываясь на графике, не критически велики)
Это квантильный график для остатков. При справедливости допущения остатки на графике должны быть близкими к наклонной прямой линии. В данном случае, остатки относительно близко прилегают к прямой (особенно в центре), а некоторые наблюдения достаточно несильно отклоняются от нее. Это говорит о корректности построенной регрессионной модели.
Если дисперсия остатков не постоянна, то облако точек будет искривленным или иметь воронкообразную форму. Здесь же график отдалённо напоминает прямую, что значит, что предположение о гомоскедастичности в целом не нарушено.
Выбросы есть, но в целом модель не объясняет их поведение и сильно подстраивается под них.
Перебирать отдельно набор переменных - долго. Воспользуемся построением полной модели.
В этой модели будет изучено влияние одновремнно всех переменных выборки на целевую.
model_3 <- lm(count ~ ., data = dplyr::select(bike_train, `holiday`:`daytime`))
summary (model_3)
##
## Call:
## lm(formula = count ~ ., data = dplyr::select(bike_train, holiday:daytime))
##
## Residuals:
## Min 1Q Median 3Q Max
## -308.1 -96.5 -30.1 55.4 686.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.91733 9.83309 0.70 0.48
## holiday -3.02912 9.71556 -0.31 0.76
## workingday 0.00856 3.47119 0.00 1.00
## weather -3.81647 2.75132 -1.39 0.17
## temp 0.48568 1.27233 0.38 0.70
## atemp 5.98645 1.17103 5.11 3.3e-07 ***
## humidity -2.19671 0.10113 -21.72 < 2e-16 ***
## windspeed 0.34126 0.21004 1.62 0.10
## month 7.31759 0.48787 15.00 < 2e-16 ***
## daytime 50.29471 1.44187 34.88 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 146 on 8698 degrees of freedom
## Multiple R-squared: 0.353, Adjusted R-squared: 0.353
## F-statistic: 528 on 9 and 8698 DF, p-value: <2e-16
Вывод: было решено сделать общую модель со всеми влиятельными переменными. Полная модель оказалась значима (p-value: <2e-16), несмотря на то, что некоторые её коэффициенты не значимы, что говорит о целесообразности её дальнейшего тестирования.
glance (model_3)
## # A tibble: 1 x 11
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
## <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int>
## 1 0.353 0.353 146. 528. 0 10 -55759. 111539. 1.12e5 1.86e8 8698
Вывод: модель опять чуть лучше справилась со своей задачей - она объяснила 35.3% всех наблюдений, и её стандартная ошибка = 146.
Построение графиков остаточной диагностики:
autoplot (model_3)
Вывод:
Рассеяние точек говорит о наличии нелинейной связи между переменными (линейная связь накладывается на пунктирную прямую в нуле), но тем не менее, коэффициенты нелинейной модели, основываясь на графике, не критически велики)
Это квантильный график для остатков. При справедливости допущения остатки на графике должны быть близкими к наклонной прямой линии. В данном случае, остатки относительно близко прилегают к прямой (особенно в центре), а некоторые наблюдения достаточно несильно отклоняются от нее. Это говорит о корректности построенной регрессионной модели.
Если дисперсия остатков не постоянна, то облако точек будет искривленным или иметь воронкообразную форму. Здесь же график отдалённо напоминает прямую, что значит, что предположение о гомоскедастичности в целом не нарушено.
Выбросы есть, но в целом модель не объясняет их поведение и сильно подстраивается под них.
Ручной перебор всех 9 переменных с целевой занял бы крайне много времени, поэтому можно оставить эту задачу компьютеру, а потом посмотреть, какие же сочетания давали наилучшие результаты. Воспользуемся информационным критерием AIc, который покажет нам, какие лучшие модели ему удалось выявитЬ.
model_empty <- lm(count ~ 1, data = dplyr::select(bike_train, `holiday`:`daytime`))
m_forward <- stepAIC(model_empty,
scope = list(lower = model_empty, upper = model_3),
direction = 'forward')
## Start: AIC=90602.3
## count ~ 1
##
## Df Sum of Sq RSS AIC
## + daytime 1 52720548 234641764 88839
## + temp 1 44730893 242631419 89131
## + atemp 1 44296295 243066017 89146
## + humidity 1 28484312 258878000 89695
## + month 1 7590617 279771696 90371
## + weather 1 5171141 282191171 90446
## + windspeed 1 2682806 284679506 90523
## <none> 287362312 90602
## + workingday 1 18893 287343420 90604
## + holiday 1 6551 287355761 90604
##
## Step: AIC=88839.3
## count ~ daytime
##
## Df Sum of Sq RSS AIC
## + temp 1 32045099 202596665 87563
## + atemp 1 31948909 202692855 87567
## + humidity 1 11973736 222668028 88385
## + month 1 8013594 226628170 88539
## + weather 1 4201890 230439874 88684
## + windspeed 1 313083 234328681 88830
## <none> 234641764 88839
## + workingday 1 27238 234614526 88840
## + holiday 1 6432 234635332 88841
##
## Step: AIC=87562.6
## count ~ daytime + temp
##
## Df Sum of Sq RSS AIC
## + humidity 1 11127784 191468881 87073
## + weather 1 3076055 199520610 87431
## + month 1 2010581 200586084 87478
## + windspeed 1 678145 201918519 87535
## + atemp 1 174827 202421838 87557
## <none> 202596665 87563
## + holiday 1 5935 202590730 87564
## + workingday 1 201 202596464 87565
##
## Step: AIC=87072.7
## count ~ daytime + temp + humidity
##
## Df Sum of Sq RSS AIC
## + month 1 5009563 186459319 86844
## + atemp 1 683656 190785226 87044
## + weather 1 161031 191307850 87067
## <none> 191468881 87073
## + windspeed 1 23090 191445791 87074
## + holiday 1 3553 191465328 87074
## + workingday 1 2411 191466470 87075
##
## Step: AIC=86843.8
## count ~ daytime + temp + humidity + month
##
## Df Sum of Sq RSS AIC
## + atemp 1 527071 185932248 86821
## + weather 1 43899 186415419 86844
## <none> 186459319 86844
## + holiday 1 5146 186454173 86846
## + windspeed 1 2790 186456528 86846
## + workingday 1 213 186459106 86846
##
## Step: AIC=86821.1
## count ~ daytime + temp + humidity + month + atemp
##
## Df Sum of Sq RSS AIC
## + windspeed 1 45077 185887171 86821
## <none> 185932248 86821
## + weather 1 29738 185902510 86822
## + holiday 1 1982 185930266 86823
## + workingday 1 23 185932225 86823
##
## Step: AIC=86821
## count ~ daytime + temp + humidity + month + atemp + windspeed
##
## Df Sum of Sq RSS AIC
## <none> 185887171 86821
## + weather 1 40979 185846192 86821
## + holiday 1 2027 185885144 86823
## + workingday 1 12 185887159 86823
Вывод: наименьшим AIC = 86821 (Информационный критерий Акаике для статистики качества приближения) обладает подобранная модель С переменными: время суток, температура,температура воздуха, влажность, месяц, температрура по ощущениям (+скорость ветра). Стоит проанализировать эту модель для проверки качества её работы.
model_4 <- lm(count ~ daytime + temp + humidity + month + atemp + windspeed, data = bike_train)
summary(model_4)
##
## Call:
## lm(formula = count ~ daytime + temp + humidity + month + atemp +
## windspeed, data = bike_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -316.4 -96.4 -29.9 55.0 689.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.6800 9.5002 0.60 0.55
## daytime 50.1097 1.4355 34.91 < 2e-16 ***
## temp 0.4384 1.2702 0.35 0.73
## humidity -2.2591 0.0907 -24.91 < 2e-16 ***
## month 7.3659 0.4865 15.14 < 2e-16 ***
## atemp 6.0360 1.1692 5.16 2.5e-07 ***
## windspeed 0.3024 0.2082 1.45 0.15
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 146 on 8701 degrees of freedom
## Multiple R-squared: 0.353, Adjusted R-squared: 0.353
## F-statistic: 792 on 6 and 8701 DF, p-value: <2e-16
Вывод: все, кроме двух коэффициентов в модели, значимые, и сама она статистически значимая p-value < 2e-16 (<0.05).
Проверка параметров модели:
glance (model_4)
## # A tibble: 1 x 11
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
## <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int>
## 1 0.353 0.353 146. 792. 0 7 -55760. 111535. 1.12e5 1.86e8 8701
Вывод: модель уже лучше справилась со своей задачей - она объяснила по-пренему объяснила 35,3% всех наблюдений, но её стандартная ошибка = 146.
Построение графиков остаточной диагностики:
autoplot (model_4)
Вывод:
Рассеяние точек говорит о наличии нелинейной связи между переменными (линейная связь накладывается на пунктирную прямую в нуле), но тем не менее, коэффициенты нелинейной модели, основываясь на графике, не критически велики)
Это квантильный график для остатков. При справедливости допущения остатки на графике должны быть близкими к наклонной прямой линии. В данном случае, остатки относительно близко прилегают к прямой (особенно в центре), а некоторые наблюдения достаточно несильно отклоняются от нее. Это говорит о корректности построенной регрессионной модели.
Если дисперсия остатков не постоянна, то облако точек будет искривленным или иметь воронкообразную форму. Здесь же график отдалённо напоминает прямую, что значит, что предположение о гомоскедастичности в целом не нарушено.
Выбросы есть, но в целом модель не объясняет их поведение и сильно подстраивается под них.
В этой модели му удалили незначимую переменную из прошлой модели и проверяем изменение работы модели после этого действия.
model_5 <- lm(count ~ daytime + humidity + month + atemp, data = bike_train)
summary(model_5)
##
## Call:
## lm(formula = count ~ daytime + humidity + month + atemp, data = bike_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -319.7 -96.9 -29.9 55.8 691.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.8944 8.2387 1.44 0.15
## daytime 50.2830 1.4306 35.15 <2e-16 ***
## humidity -2.2988 0.0872 -26.37 <2e-16 ***
## month 7.3071 0.4848 15.07 <2e-16 ***
## atemp 6.4143 0.1941 33.05 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 146 on 8703 degrees of freedom
## Multiple R-squared: 0.353, Adjusted R-squared: 0.353
## F-statistic: 1.19e+03 on 4 and 8703 DF, p-value: <2e-16
Вывод: все коэффициенты в модели значимые, и сама она статистически значимая p-value < 2e-16 (<0.05).
Проверка параметров модели:
glance (model_5)
## # A tibble: 1 x 11
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
## <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int>
## 1 0.353 0.353 146. 1187. 0 5 -55761. 111534. 1.12e5 1.86e8 8703
Вывод: модель также справилась со своей задачей - она объяснила по-пренему объяснила 35,3% всех наблюдений, но её стандартная ошибка = 146.
Построение графиков остаточной диагностики:
autoplot (model_5)
Вывод:
Рассеяние точек говорит о наличии нелинейной связи между переменными (линейная связь накладывается на пунктирную прямую в нуле), но тем не менее, коэффициенты нелинейной модели, основываясь на графике, не критически велики)
Это квантильный график для остатков. При справедливости допущения остатки на графике должны быть близкими к наклонной прямой линии. В данном случае, остатки относительно близко прилегают к прямой (особенно в центре), а некоторые наблюдения достаточно несильно отклоняются от нее. Это говорит о корректности построенной регрессионной модели.
Если дисперсия остатков не постоянна, то облако точек будет искривленным или иметь воронкообразную форму. Здесь же график отдалённо напоминает прямую, что значит, что предположение о гомоскедастичности в целом не нарушено.
Выбросы есть, и они хуже объясняются, чем в предыдущий модели.
Теперь отойдём от применения простых регрессионных моделей и рассмотрим более сложные. Для визуализации нескольких регресионных прямых, преобразуем рассамтриваемые переменные в отдельных данных в качественные:
bike_train1 <- bike_train
bike_train1$daytime <- as.character(bike_train1$daytime)
bike_train1 <- bike_train1
bike_train1$weather <- as.character(bike_train1$weather)
bike_train1 <- bike_train1
bike_train1$month <- as.character(bike_train1$month)
bike_train1 <- bike_train1
bike_train1$workingday <- as.character(bike_train1$workingday)
Теперь построим визуализации регресионных прямых:
dots1 <- bike_train1 %>%
ggplot(aes(count, temp, color = workingday)) +
geom_point(alpha = 0.1) +
geom_smooth(method = 'lm', se = F)+
labs(title="1", colour = "Рабочие дни", x = "Число актов аренды", y = "Температура воздуха")
dots2 <-bike_train1 %>%
ggplot(aes(count, temp, color = daytime)) +
geom_point(alpha = 0.1) +
geom_smooth(method = 'lm', se = F)+
labs(title="2",colour = "Время суток", x = "Число актов аренды", y = "Температура воздуха")
dots3 <-bike_train1 %>%
ggplot(aes(count, temp, color = month)) +
geom_point(alpha = 0.1) +
geom_smooth(method = 'lm', se = F)+
labs(title="3", colour = "Месяц", x = "Число актов аренды", y = "Температура воздуха")
dots4 <-bike_train1 %>%
ggplot(aes(count, temp, color = weather)) +
geom_point(alpha = 0.1) +
geom_smooth(method = 'lm', se = F)+
labs(title="4",colour = "Состояние погоды", x = "Число актов аренды", y = "Температура воздуха")
grid.arrange(dots1, dots2, dots3, dots4, ncol=2)
Вывод:
1 график: Какой-либо зависимости между переменными незамечено. Можно сказать, что происходит много числов актов аренды при высокой температуре воздуха в выходные данные, но данная закономерность выражена слаба. Всё же никаких взаимосвязей между этими переменными обнаружено не было;
2 график: Можно заявить, что независимо от температуры воздуха, в ночное время суток берётся в аренду крайне низкое число велосипедов. В остальные времена суток - распределение актов аренды равномерно и не зависит от температуры воздуха;
3 график: Можно заметить, что в зимние времена года, когда присутствовует низкая температура воздуха, число актов аренды минимально. В остальных месяцах число актов аренды распределено равномерно без подчинения к температуре воздуха;
4 график: Сразу бросается в глаза тот факт, что большинство актов аренды происходят во время погоды, при которой ясно, малооблачно или немного облачно. Крайне мало велосипедов арендавалось во время ливня, града, грома, тумана и снега. Другими словами, состояние погоды влияет на число актов аренды по ухудшению ей условий (1 тип -наибольшее число, 4 тип - наименьшее число);
model_6<- lm(count ~ temp + daytime + humidity + daytime*workingday + month + weather*daytime, data = bike_train)
summary(model_6)
##
## Call:
## lm(formula = count ~ temp + daytime + humidity + daytime * workingday +
## month + weather * daytime, data = bike_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -324.6 -97.2 -30.3 55.0 693.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.5730 12.3121 0.62 0.53851
## temp 6.8673 0.2114 32.48 < 2e-16 ***
## daytime 55.8322 3.8841 14.37 < 2e-16 ***
## humidity -2.1476 0.0967 -22.20 < 2e-16 ***
## workingday -27.8581 7.8322 -3.56 0.00038 ***
## month 7.3488 0.4860 15.12 < 2e-16 ***
## weather 17.6532 5.7794 3.05 0.00226 **
## daytime:workingday 11.2582 2.9147 3.86 0.00011 ***
## daytime:weather -8.9447 2.1408 -4.18 3e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 146 on 8699 degrees of freedom
## Multiple R-squared: 0.354, Adjusted R-squared: 0.353
## F-statistic: 595 on 8 and 8699 DF, p-value: <2e-16
Вывод: все коэффициенты в модели значимые, и сама она статистически значимая p-value < 2e-16 (<0.05).
Проверка параметров работы моделей:
glance(model_6)
## # A tibble: 1 x 11
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
## <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int>
## 1 0.354 0.353 146. 595. 0 9 -55756. 111532. 1.12e5 1.86e8 8699
Вывод: модель объяснила 35.4% наблюдений, а стандартное отклонение = 146. Пока что эта модель с взаимосвязями между переменными показала наилучшие результаты.
Диагностические графики:
autoplot (model_6)
Вывод:
Рассеяние точек говорит о наличии нелинейной связи между переменными (линейная связь накладывается на пунктирную прямую в нуле), но тем не менее, коэффициенты нелинейной модели, основываясь на графике, не критически велики).
Это квантильный график для остатков. При справедливости допущения остатки на графике должны быть близкими к наклонной прямой линии. В данном случае, остатки относительно близко прилегают к прямой (особенно в центре), а некоторые наблюдения достаточно несильно отклоняются от нее. Это говорит о корректности построенной регрессионной модели.
Если дисперсия остатков не постоянна, то облако точек будет искривленным или иметь воронкообразную форму. Здесь же график отдалённо напоминает прямую, что значит, что предположение о гомоскедастичности в целом не нарушено.
Левередж = 3*10/8708 = 0.0035 Выбросы есть, но в целом модель несильно подстраивается под них. Достаточно много наблюдений выходят за лопустимые границы, т.е. являются весомыми при построении модели. Возможно, стоит произвести анализ выбросов.
model_6a <- model_6 %>%
fortify() %>% # сохранили остатки, предсказания и показатели влиятельности в набор
merge(., bike_train) # добавили исходные данные (merge удаляет дублирующиеся при слиянии столбцы)
ggplot(model_6a, aes(.fitted, count)) +
geom_point(aes(size = .cooksd, colour = .hat)) +
geom_line(aes(.fitted, .fitted),
colour = "red", linetype = "dashed") +
scale_colour_gradient(low = "black", high = "red") +
ggrepel::geom_label_repel(
data = filter(model_6a, .cooksd > 4/nrow(bike_train)),
mapping = aes(label = "danger"), alpha = 0.25) +
labs(title = "Влиятельные наблюдения",
x = "Прогноз count", y = "Факт count",
colour = "Левередж", size = "Расстояние Кука")
Влиятельные наблюдения помечены словом “danger”, что говорит об опасности использования данных наблюдений в моделе, которые негавтивно на нее влияют. Уберем выбросы и пересчитаем лучшие прогнозы и ошибки моделей:
model_6a$.cooksd <- as.numeric(model_6a$.cooksd)
model_6a$.cooksd[is.na(model_6a$.cooksd)] <- mean(model_6a$.cooksd, na.rm = T)
bike_train_clean <- model_6a %>%
filter(.cooksd < 4/nrow(model_6a))
Определение степени целевой переменной для устранения гетеродоксичности:
t3<- car::powerTransform(count ~ temp + daytime + humidity + daytime*workingday + month + weather*daytime, data = bike_train_clean)
coef(t3)
## Y1
## 0.308175
Построение модели:
model_6b<- lm(count^0.308 ~ temp + daytime + I(daytime^2)+ workingday + I(month^2)+ humidity + daytime*workingday + month + weather*daytime, data = bike_train_clean)
summary(model_6b)
##
## Call:
## lm(formula = count^0.308 ~ temp + daytime + I(daytime^2) + workingday +
## I(month^2) + humidity + daytime * workingday + month + weather *
## daytime, data = bike_train_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.254 -0.669 -0.051 0.603 3.112
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.540505 0.102399 -15.04 < 2e-16 ***
## temp 0.035696 0.002358 15.14 < 2e-16 ***
## daytime 4.206833 0.061073 68.88 < 2e-16 ***
## I(daytime^2) -0.695317 0.010719 -64.87 < 2e-16 ***
## workingday -0.438254 0.048477 -9.04 < 2e-16 ***
## I(month^2) -0.011797 0.001643 -7.18 7.7e-13 ***
## humidity -0.008818 0.000666 -13.23 < 2e-16 ***
## month 0.226562 0.022873 9.91 < 2e-16 ***
## weather 0.029493 0.036046 0.82 0.41
## daytime:workingday 0.136210 0.018206 7.48 8.1e-14 ***
## daytime:weather -0.079945 0.013671 -5.85 5.2e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.907 on 8443 degrees of freedom
## Multiple R-squared: 0.684, Adjusted R-squared: 0.684
## F-statistic: 1.83e+03 on 10 and 8443 DF, p-value: <2e-16
Вывод: все коэффициенты в модели, кроме одного, значимые, и сама она статистически значимая p-value < 2e-16 (<0.05).
Диагностические графики:
autoplot(model_6b)
Вывод: практически все диагностические графики демонстрируют превосходную работу модели. Самое главное, что нам удалось добиться практически полной горизонтальной прямой на 3-ем графике.
Поменяем в моделе без выбросов переменную temp на переменную atemp, которая обладала большим коэффициентом корреляции при разведоточном анализе:
t4<- car::powerTransform(count ~ atemp + daytime + humidity + daytime*workingday + month + weather*daytime, data = bike_train_clean)
coef(t4)
## Y1
## 0.308446
Построение модели:
model_6b1<- lm(count^0.308 ~ atemp + daytime+I(daytime^2) + humidity + daytime*workingday + month + weather*daytime, data = bike_train_clean)
summary(model_6b1)
##
## Call:
## lm(formula = count^0.308 ~ atemp + daytime + I(daytime^2) + humidity +
## daytime * workingday + month + weather * daytime, data = bike_train_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.214 -0.667 -0.055 0.606 3.179
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.373797 0.098361 -13.97 < 2e-16 ***
## atemp 0.045904 0.001235 37.17 < 2e-16 ***
## daytime 4.125764 0.060305 68.41 < 2e-16 ***
## I(daytime^2) -0.683799 0.010628 -64.34 < 2e-16 ***
## humidity -0.008710 0.000662 -13.16 < 2e-16 ***
## workingday -0.436680 0.048559 -8.99 < 2e-16 ***
## month 0.063237 0.003063 20.65 < 2e-16 ***
## weather 0.009930 0.036011 0.28 0.78
## daytime:workingday 0.136447 0.018239 7.48 8.1e-14 ***
## daytime:weather -0.072265 0.013682 -5.28 1.3e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.909 on 8444 degrees of freedom
## Multiple R-squared: 0.683, Adjusted R-squared: 0.682
## F-statistic: 2.02e+03 on 9 and 8444 DF, p-value: <2e-16
Вывод: все коэффициенты в модели, кроме одного, значимые, и сама она статистически значимая p-value < 2e-16 (<0.05).
Определеим степень, в которую необходимо возвести целевую переменную, чтобы сделать нашу лучшую модель гомодоксичной:
t2<- car::powerTransform(count ~ atemp + daytime + humidity + daytime*workingday + month*daytime + weather*daytime + windspeed + humidity*month, data = bike_train_clean)
coef(t2)
## Y1
## 0.314036
Построение модели:
model_6b2m<- lm(count^0.314 ~ atemp + daytime + I(daytime^2) + humidity + daytime*workingday + month*daytime + weather*daytime + humidity*month + humidity*daytime, data = bike_train_clean)
summary(model_6b2m)
##
## Call:
## lm(formula = count^0.314 ~ atemp + daytime + I(daytime^2) + humidity +
## daytime * workingday + month * daytime + weather * daytime +
## humidity * month + humidity * daytime, data = bike_train_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.436 -0.698 -0.062 0.630 3.427
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.889901 0.146525 -12.90 < 2e-16 ***
## atemp 0.047121 0.001298 36.30 < 2e-16 ***
## daytime 4.349183 0.072610 59.90 < 2e-16 ***
## I(daytime^2) -0.717278 0.011201 -64.04 < 2e-16 ***
## humidity 0.001762 0.001893 0.93 0.3521
## workingday -0.456241 0.050677 -9.00 < 2e-16 ***
## month 0.087614 0.014039 6.24 4.6e-10 ***
## weather -0.052298 0.040383 -1.30 0.1953
## daytime:workingday 0.142583 0.019033 7.49 7.5e-14 ***
## daytime:month 0.013635 0.002766 4.93 8.4e-07 ***
## daytime:weather -0.048311 0.015746 -3.07 0.0022 **
## humidity:month -0.000858 0.000166 -5.16 2.5e-07 ***
## daytime:humidity -0.002319 0.000591 -3.92 8.8e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.949 on 8441 degrees of freedom
## Multiple R-squared: 0.685, Adjusted R-squared: 0.684
## F-statistic: 1.53e+03 on 12 and 8441 DF, p-value: <2e-16
Вывод: не все коэффициенты в модели значимые, но сама она статистически значимая p-value < 2e-16 (<0.05).
autoplot(model_6b2m)
Вывод: практически все диагностические графики демонстрируют превосходную работу модели. Самое главное, что нам удалось добиться практически полной горизонтальной прямой на 3-ем графике.
Подключение необходимых библиотек для работы с деревьями решений:
library(scales) # Процентный формат для осей ggplot
##
## Attaching package: 'scales'
## The following object is masked from 'package:memisc':
##
## percent
## The following object is masked from 'package:purrr':
##
## discard
## The following object is masked from 'package:readr':
##
## col_factor
library(rpart) # Деревья в R: recursive partitioning
library(party) # Деревья в R: conditional inference
## Loading required package: grid
## Loading required package: mvtnorm
## Loading required package: modeltools
## Loading required package: stats4
##
## Attaching package: 'modeltools'
## The following objects are masked from 'package:memisc':
##
## Lapply, relabel
## The following object is masked from 'package:car':
##
## Predict
## Loading required package: strucchange
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sandwich
##
## Attaching package: 'strucchange'
## The following object is masked from 'package:stringr':
##
## boundary
library(rpart.plot) # Визуализация деревьев
library(partykit) # Визуализация деревьев
## Loading required package: libcoin
##
## Attaching package: 'partykit'
## The following objects are masked from 'package:party':
##
## cforest, ctree, ctree_control, edge_simple, mob, mob_control, node_barplot,
## node_bivplot, node_boxplot, node_inner, node_surv, node_terminal, varimp
library(C50) # Деревья в R: C5.0
library(mlr)
## Loading required package: ParamHelpers
## 'mlr' is in maintenance mode since July 2019. Future development efforts will go into its
## successor 'mlr3' (<https://mlr3.mlr-org.com>).
##
## Attaching package: 'mlr'
## The following object is masked from 'package:modelr':
##
## resample
Построение автоматически подобронного дерева решений на обучающей выборке:
rtree <- rpart(formula = count~., data = bike_train)
summary(rtree)
## Call:
## rpart(formula = count ~ ., data = bike_train)
## n= 8708
##
## CP nsplit rel error xerror xstd
## 1 0.3108764 0 1.000000 1.000241 0.0195569
## 2 0.0907987 1 0.689124 0.689303 0.0155017
## 3 0.0263655 2 0.598325 0.603862 0.0137338
## 4 0.0213274 3 0.571959 0.574431 0.0131302
## 5 0.0114109 4 0.550632 0.561840 0.0129384
## 6 0.0100000 5 0.539221 0.556046 0.0127668
##
## Variable importance
## daytime temp atemp month humidity weather
## 53 15 15 10 5 1
##
## Node number 1: 8708 observations, complexity param=0.310876
## mean=191.136, MSE=32999.8
## left son=2 (2523 obs) right son=3 (6185 obs)
## Primary splits:
## daytime < 1.5 to the left, improve=0.3108760, (0 missing)
## atemp < 29.925 to the left, improve=0.1320830, (0 missing)
## temp < 22.55 to the left, improve=0.1044840, (0 missing)
## humidity < 66.5 to the right, improve=0.0716794, (0 missing)
## month < 3.5 to the left, improve=0.0548867, (0 missing)
## Surrogate splits:
## temp < 4.51 to the left, agree=0.712, adj=0.006, (0 split)
## atemp < 3.41 to the left, agree=0.711, adj=0.004, (0 split)
##
## Node number 2: 2523 observations
## mean=32.5513, MSE=1584.38
##
## Node number 3: 6185 observations, complexity param=0.0907987
## mean=255.826, MSE=31371.2
## left son=6 (2717 obs) right son=7 (3468 obs)
## Primary splits:
## temp < 19.27 to the left, improve=0.1344740, (0 missing)
## atemp < 23.105 to the left, improve=0.1334120, (0 missing)
## month < 3.5 to the left, improve=0.1077000, (0 missing)
## humidity < 74.5 to the right, improve=0.0399052, (0 missing)
## weather < 2.5 to the right, improve=0.0275818, (0 missing)
## Surrogate splits:
## atemp < 23.105 to the left, agree=0.998, adj=0.996, (0 split)
## month < 3.5 to the left, agree=0.751, adj=0.434, (0 split)
## humidity < 85.5 to the right, agree=0.578, adj=0.040, (0 split)
## windspeed < 21.0012 to the right, agree=0.567, adj=0.014, (0 split)
##
## Node number 6: 2717 observations, complexity param=0.0213274
## mean=182.446, MSE=19869.5
## left son=12 (1665 obs) right son=13 (1052 obs)
## Primary splits:
## month < 9.5 to the left, improve=0.1135250, (0 missing)
## atemp < 14.015 to the left, improve=0.0824801, (0 missing)
## temp < 12.71 to the left, improve=0.0817189, (0 missing)
## daytime < 3.5 to the right, improve=0.0431468, (0 missing)
## weather < 2.5 to the right, improve=0.0259138, (0 missing)
## Surrogate splits:
## humidity < 93.5 to the left, agree=0.617, adj=0.011, (0 split)
## atemp < 17.045 to the left, agree=0.613, adj=0.001, (0 split)
##
## Node number 7: 3468 observations, complexity param=0.0263655
## mean=313.316, MSE=32858.5
## left son=14 (875 obs) right son=15 (2593 obs)
## Primary splits:
## humidity < 72.5 to the right, improve=0.0664875, (0 missing)
## atemp < 29.925 to the left, improve=0.0459233, (0 missing)
## weather < 2.5 to the right, improve=0.0327580, (0 missing)
## temp < 24.19 to the left, improve=0.0226283, (0 missing)
## daytime < 2.5 to the left, improve=0.0173916, (0 missing)
## Surrogate splits:
## weather < 2.5 to the right, agree=0.795, adj=0.187, (0 split)
##
## Node number 12: 1665 observations
## mean=144.694, MSE=14110.9
##
## Node number 13: 1052 observations
## mean=242.196, MSE=23157.8
##
## Node number 14: 875 observations
## mean=232.854, MSE=24604.7
##
## Node number 15: 2593 observations, complexity param=0.0114109
## mean=340.467, MSE=32721.8
## left son=30 (749 obs) right son=31 (1844 obs)
## Primary splits:
## daytime < 2.5 to the left, improve=0.0386463, (0 missing)
## humidity < 44.5 to the right, improve=0.0263598, (0 missing)
## temp < 23.37 to the left, improve=0.0228858, (0 missing)
## atemp < 26.895 to the left, improve=0.0224512, (0 missing)
## workingday < 0.5 to the right, improve=0.0157984, (0 missing)
##
## Node number 30: 749 observations
## mean=284.67, MSE=21002.9
##
## Node number 31: 1844 observations
## mean=363.131, MSE=35703.6
rpart.plot(rtree)
Вывод: верхнее число над процентом (191, 256 и т.д.) - это среднее число арендованных велосипедов, исходя из условий отбора. Данная модель, по сути, способна прогнозировать только 6 значений count для каждого дня. Модель самостоятельно отобрала, по её мнению, самые важные/значимые переменные - daytime, temp, humidity, month. Примерно такие же переменные были отобраны при ручном моделировании. Значит, данной моделе можно доверять. Попробуем увеличить число ветвей - это позволит повысить число прогнозируемых значений, что в теории должно повысить точность. Тем не менее, позже посмтроим работу этой модели на обучающей выборке.
Вычисление среднего значения арендованных велосипедов за весь обучающий период:
mean(bike_train$count)
## [1] 191.136
Создание дополнительного столбца с разделение в числах актов аренды велосипедов по признаку “больше или меньшее, чем среднее”:
bike_train_a <- bike_train %>%
mutate( count1 = ifelse(count>=191, "much", "less"))
График рассеяния актов аренды велосипедов в разрезе времени суток, температуры и с разделением по признаку “больше или меньшее, чем среднее”:
cmap <- c("much" = "green", "less" = "red")
ggplot(bike_train_a, aes(x = daytime, y = temp)) +
geom_jitter(aes(colour = count1), alpha = 0.5) +
labs(title = "Сравнение числа актов аренды велосипедов и его среднего значения", colour = "Число аренды в день", x = "Время суток", y = "t воздуха") +
scale_x_continuous(breaks = 1:10) +
scale_colour_manual(values = cmap)
Вывод: можно заметить, что ночью практически каждое наблюдение содержит число актов аренды меньшее, чем среднее значение за весь обучащий период, а также низкая температура воздуха, независимо от времени суток негативно влияет на желание людей брать велосипеды и они берут их реже обычного.
Построение модифицированного дерева решений:
rtree_overfit <- rpart(count ~ .,
data = bike_train,
control =
rpart.control(
minsplit = 2, # минимальное число примеров для разбиения
cp = 0.00001, # минимальное относительное улучшение для разбиения
maxdepth = 20 # максимальная глубина для листа
))
rpart.plot(rtree_overfit)
Вывод: модель получилась чересчур сложной, она нуждается в упращении.
Посмторим на зависимость относительной ошибки от cp:
plotcp(rtree_overfit)
Вывод: можно увидеть, как сложность модели влияет на ее точность. Попробуем установить параметр сложности (ср) равный 0.0003.
Древо решений с оптимальной cp:
rtree_overfit1 <- prune(rtree_overfit, cp = 0.0003)
rpart.plot(rtree_overfit1)
Вывод: модель стала проще, но тем не меннее, по-прежнему сложна. Также её прогноза на обучающей выборке
Посмотрим на полную дерево решений по обучающей выборке:
rctree <- ctree(count ~ ., data = bike_train)
plot(rctree)
Вывод: она оказалась так же сложна. На всякий случай проверим точность её прогноза на обучающей выборке.
Подключение необходимых библиотек для работы с автоматическими моделями:
library(tidyverse) # трансформация и визуализация данных
library(mlr) # фреймворк для машинного обучения
library(stringr) # работа со строками
library(forcats) # работа с факторами
library(parallelMap) # распараллеливание задач
library(rpart.plot) # визуализация деревьев rpart
Созданим задачу для автоматической модели: целевая переменная count, данные bike:
bike_task <- makeRegrTask(id = "Аренда велосипедов",
data = bike, target = "count",
fixup.data = "warn")
print(bike_task)
## Supervised task: Аренда велосипедов
## Type: regr
## Target: count
## Observations: 10886
## Features:
## numerics factors ordered functionals
## 9 0 0 0
## Missings: FALSE
## Has weights: FALSE
## Has blocking: FALSE
## Has coordinates: FALSE
Посмотрим, какие модели существуют для автоматического прогнозирования на основе регресии:
listLearners(obj = "regr") %>%
dplyr::select(name, class) %>%
mutate(name = substr(name, 1, 60)) %>%
head(10) %>%
arrange(name)
## name class
## 1 Bayesian Additive Regression Trees regr.bartMachine
## 2 Bayesian CART regr.bcart
## 3 Bayesian Gaussian Process regr.bgp
## 4 Bayesian Gaussian Process with jumps to the Limiting Linear regr.bgpllm
## 5 Bayesian Linear Model regr.blm
## 6 Bayesian regularization for feed-forward neural networks regr.brnn
## 7 Bayesian Treed Gaussian Process regr.btgp
## 8 Bayesian Treed Gaussian Process with jumps to the Limiting L regr.btgpllm
## 9 Bayesian Treed Linear Model regr.btlm
## 10 Gradient Boosting regr.bst
Поиск точного колиества таких моделей
listLearners(obj = "regr") %>% count()
## # A tibble: 1 x 1
## n
## <int>
## 1 59
Вывод: их количество - 59.
Прогноз по самым популярным и точным моделям: линейная регрессия:
lrn_lm <- makeLearner("regr.lm", id = "Линейная регрессия")
lrn_lm
## Learner Линейная регрессия from package stats
## Type: regr
## Name: Simple Linear Regression; Short name: lm
## Class: regr.lm
## Properties: numerics,factors,se,weights
## Predict-Type: response
## Hyperparameters:
Взаимосвязь коэффициентов модели и её проверка:
m_lm <- train(lrn_lm, bike_task)
m_lm %>%
getLearnerModel() %>%
summary()
##
## Call:
## stats::lm(formula = f, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -314.2 -96.5 -30.7 55.9 685.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.0220 8.7741 1.26 0.209
## holiday -4.6855 8.6788 -0.54 0.589
## workingday -0.0405 3.1053 -0.01 0.990
## weather -3.4616 2.4667 -1.40 0.161
## temp 2.3462 1.0698 2.19 0.028 *
## atemp 4.2313 0.9846 4.30 1.7e-05 ***
## humidity -2.1815 0.0902 -24.19 < 2e-16 ***
## windspeed 0.3420 0.1876 1.82 0.068 .
## month 7.3987 0.4343 17.03 < 2e-16 ***
## daytime 49.6924 1.2903 38.51 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 146 on 10876 degrees of freedom
## Multiple R-squared: 0.352, Adjusted R-squared: 0.351
## F-statistic: 656 on 9 and 10876 DF, p-value: <2e-16
Вывод: некоторые переменные не значимы, но в целом, модель можно использовать, так как сама она статистичсеки значима?: p-value<0.05.
Разделение данных на обучающий и тестовый период для работ с автоматическими моделями:
# Получаем номера строк для обучающей и тестовой выборки
set.seed(123) # для воспроизводимости
bike_n <- getTaskSize(bike_task)
bike_train_task <- sample(bike_n, size = bike_n * 0.7)
bike_test_task <- base::setdiff(1:bike_n, bike_train_task)
# Строим модель только на обучающей выборке
bike_mod_lm <- train(lrn_lm, task = bike_task,
subset = bike_train_task)
Прогноз с помощью неё на обучающем периоде:
bike_pred_lm <- predict(bike_mod_lm,
task = bike_task,
subset = bike_train_task)
bike_pred_lm
## Prediction: 7620 observations
## predict.type: response
## threshold:
## time: 0.00
## id truth response
## 2463 2463 207 298.956
## 2511 2511 420 293.396
## 10419 10419 242 178.309
## 8718 8718 526 402.898
## 2986 2986 313 264.311
## 1842 1842 154 277.009
## ... (#rows: 7620, #cols: 3)
Расчет ошибок модели:
reg_ms <- list(rmse, mape, mae)
performance(bike_pred_lm,
measures = reg_ms) %>%
round(2)
## rmse mape mae
## 144.95 3.16 107.44
Вывод:точность модели плохая, поэтому будут использованы альтернативные автоматические модели для их последующего включение в отбор.
lrn_lm1 <- makeLearner("regr.randomForest")
lrn_lm1
## Learner regr.randomForest from package randomForest
## Type: regr
## Name: Random Forest; Short name: rf
## Class: regr.randomForest
## Properties: numerics,factors,ordered,se,oobpreds,featimp
## Predict-Type: response
## Hyperparameters:
Присвоение модели выборки:
m_lm1 <- train(lrn_lm1, bike_task)
Просмотр показателей модели (переменные), которые сложно интерпретировать:
m_lm1 %>%
getLearnerModel() %>%
summary()
## Length Class Mode
## call 4 -none- call
## type 1 -none- character
## predicted 10886 -none- numeric
## mse 500 -none- numeric
## rsq 500 -none- numeric
## oob.times 10886 -none- numeric
## importance 9 -none- numeric
## importanceSD 0 -none- NULL
## localImportance 0 -none- NULL
## proximity 0 -none- NULL
## ntree 1 -none- numeric
## mtry 1 -none- numeric
## forest 11 -none- list
## coefs 0 -none- NULL
## y 10886 -none- numeric
## test 0 -none- NULL
## inbag 0 -none- NULL
Разделение данных на обучающую и тестовую выборки и проноз:
# Получаем номера строк для обучающей и тестовой выборки
set.seed(123) # для воспроизводимости
bike_n1 <- getTaskSize(bike_task)
bike_train_task1 <- sample(bike_n1, size = bike_n1 * 0.7)
bike_test_task1 <- base::setdiff(1:bike_n1, bike_train_task1)
# Строим модель только на обучающей выборке
bike_mod_lm1 <- train(lrn_lm1, task = bike_task,
subset = bike_train_task1)
bike_pred_lm1 <- predict(bike_mod_lm1,
task = bike_task,
subset = bike_train_task1)
bike_pred_lm1
## Prediction: 7620 observations
## predict.type: response
## threshold:
## time: 0.53
## id truth response
## 2463 2463 207 202.916
## 2511 2511 420 392.187
## 10419 10419 242 262.529
## 8718 8718 526 472.980
## 2986 2986 313 333.613
## 1842 1842 154 201.492
## ... (#rows: 7620, #cols: 3)
Расчет ошибок модели на обучающем периоде:
reg_ms1 <- list(rmse, mape, mae)
performance(bike_pred_lm1,
measures = reg_ms) %>%
round(2)
## rmse mape mae
## 66.60 0.80 44.57
Вывод: Малые ошибки модели связаны с тем, что randomForest в результате обучения испльзует не одно, а несколько деревьев решений, которые по отдельности показывают не лучшие результаты. Точность получилась неплохой не из-за того, что данный метод для оценки хороший, а за счет чересчур множественных ответвлений и условий при построении дерева. Возможно, что тестовый период окажется не таким точным.
Модель уже лучше, по сравнению с lm, но скорее всего, уступает ручным моделям. Тем не менее, она поучаствует в конкурсе.
*Далее, мы решили посчитать проноз по модели “нечетких множеств”. В течение нашего обучения мы часто сталкивались с данным подходом, который помогал решать множесто логистических задач, связанных с неопределенностью. Здесь не будут представлены chunks с этой моделью, т.к. ее выполнение требует очень долгого времени и в целом нагружает компьютер. Однако, ошибки модели будут представлены: rmse: 132.58 mape: 1.87 mae: 94.95 Результаты ошибок не самые лучшие, значит, модель не будет участвовать в борьбе за лучшую
MLR содержит 59 моделей регресии, считать каждую вручную - нецелесообразно. Попробуем отобрать модели, наиболее подходящие к исследуемым данным. Исходя из характеристик, которыми могут быть описаны модели, выделим 2 значимые: 1) SE - позволяет рассчитывать и прогнозировать стандартные ошибки модели. это позволит выбирать лучшую модель. 2) Numeric - позволяет прогнозировать числовые переменные (count - целевая числовая переменная) В остальном данные не имеют каких-либо особых харакетристик, которые позволили бы отобрать подходящие модели.
listLearners(obj = "regr", properties = c("se", "numerics")) %>%
dplyr::select(name, class) %>%
mutate(name = substr(name, 1, 60)) %>%
head(50) %>%
arrange(name)
## name class
## 1 Bayesian CART regr.bcart
## 2 Bayesian Gaussian Process regr.bgp
## 3 Bayesian Gaussian Process with jumps to the Limiting Linear regr.bgpllm
## 4 Bayesian Linear Model regr.blm
## 5 Bayesian Treed Gaussian Process regr.btgp
## 6 Bayesian Treed Gaussian Process with jumps to the Limiting L regr.btgpllm
## 7 Bayesian Treed Linear Model regr.btlm
## 8 Gaussian Process regr.GPfit
## 9 Gaussian Processes regr.gausspr
## 10 Generalized Linear Regression regr.glm
## 11 Kriging regr.km
## 12 Local Approximate Gaussian Process regr.laGP
## 13 Random Forest regr.randomForest
## 14 Random Forests regr.ranger
## 15 Regression Splines regr.crs
## 16 Simple Linear Regression regr.lm
Вывод: итого было отобрано 16 моделей, соответствующие заданным параметрам.
Подключение необходимых библиотек для работы с отобранными моделями:
library(tgp)
library(GPfit)
library(DiceKriging)
library(laGP)
library(ranger)
library(crs)
## Registered S3 method overwritten by 'crs':
## method from
## predict.gsl.bs np
## Categorical Regression Splines (version 0.15-31)
## [vignette("crs_faq") provides answers to frequently asked questions]
Вывод: lm и randomforest уже расчитаны. Уберем их.
Создание нескольких моделей:
turnover_learners <-
makeLearners(cls = c("bcart", "bgp", "bgpllm", "blm", "btgp", "btgpllm", "btlm", "GPfit","gausspr", "glm", "km", "laGP", "ranger", "crs" ),
ids = c("bcart", "bgp", "bgpllm", "blm", "btgp", "btgpllm", "btlm", "GPfit","gausspr", "glm", "km", "laGP", "ranger", "crs"),
type = "regr", predict.type = "se")
Создаем задачу - обозначаем целевую переменную:
turnover_task <- makeRegrTask(id = "Аренда велосипедов",
data = bike, target = "count",
fixup.data = "warn")
Облегчение процесса расчета прогнозов:
set.seed(123, "L'Ecuyer") # При использовании параллельных вычислений необходимо выбирать алгоритм генерации СЧ "L'Ecuyer"
# Определение количества доступных ядер процессора
num_cores <- parallel::detectCores()
# Запуск параллельных вычислений
parallelStartSocket(num_cores) # для Windows
## Starting parallelization in mode=socket with cpus=8.
Вывод прогнозирование (Данный chunk не выполняется до конца - видимо требует слишком много оперативной памяти, компьютеры не справляются. При уменьшении числа моделей в makelearners ситуация не меняется. Было принято решение - рассчитать модели вручную по отдельности. Ошибки моделей будут представлены ниже.)
turnover_bench <- turnover_learners %>% benchmark(tasks = turnover_task, show.info = FALSE)
“bcart”, “bgp”, “bgpllm”, “blm”, “btgp”, “btgpllm”, “btlm”, “GPfit”,“gausspr”, “glm”, “km”, “laGP”, “ranger”, “crs”
Постановка задачи моделям:
lrn_bcart <- makeLearner("regr.bcart")
lrn_bcart
## Learner regr.bcart from package tgp
## Type: regr
## Name: Bayesian CART; Short name: bcart
## Class: regr.bcart
## Properties: numerics,se,factors
## Predict-Type: response
## Hyperparameters:
Присвоение моделям выборки:
m_bcart1 <- train(lrn_lm1, bike_task)
Показатели модели (переменные), котоыре сложно интерпретировать:
m_bcart1 %>%
getLearnerModel() %>%
summary()
## Length Class Mode
## call 4 -none- call
## type 1 -none- character
## predicted 10886 -none- numeric
## mse 500 -none- numeric
## rsq 500 -none- numeric
## oob.times 10886 -none- numeric
## importance 9 -none- numeric
## importanceSD 0 -none- NULL
## localImportance 0 -none- NULL
## proximity 0 -none- NULL
## ntree 1 -none- numeric
## mtry 1 -none- numeric
## forest 11 -none- list
## coefs 0 -none- NULL
## y 10886 -none- numeric
## test 0 -none- NULL
## inbag 0 -none- NULL
Разделение данных на обучающую и тестовую выборки и прогноз:
# Получаем номера строк для обучающей и тестовой выборки
set.seed(123) # для воспроизводимости
bike_n1 <- getTaskSize(bike_task)
bike_train_task1 <- sample(bike_n1, size = bike_n1 * 0.7)
bike_test_task1 <- base::setdiff(1:bike_n1, bike_train_task1)
# Строим модель только на обучающей
bike_mod_bcart1 <- train(lrn_bcart, task = bike_task,
subset = bike_train_task1)
##
## burn in:
## **GROW** @depth 0: [5,0.499943], n=(3619,4001)
## **GROW** @depth 1: [5,0.706953], n=(2530,1340)
## **GROW** @depth 2: [5,0.793103], n=(1576,434)
## **GROW** @depth 2: [6,0.43], n=(299,1965)
## **GROW** @depth 3: [3,0], n=(1172,404)
## **GROW** @depth 2: [2,0], n=(135,299)
## **GROW** @depth 2: [6,0.5], n=(983,1889)
## **PRUNE** @depth 2: [2,0]
## **GROW** @depth 3: [9,0], n=(78,326)
## **GROW** @depth 2: [9,0.333333], n=(1208,448)
## **GROW** @depth 3: [8,0.454545], n=(726,482)
## **GROW** @depth 3: [2,0], n=(174,406)
## **GROW** @depth 5: [4,0.489796], n=(580,1041)
## **GROW** @depth 2: [2,0], n=(428,850)
## **GROW** @depth 5: [6,0.56], n=(180,1126)
## **GROW** @depth 3: [6,0.49], n=(176,252)
## **PRUNE** @depth 3: [6,0.53]
## **GROW** @depth 3: [4,0.285714], n=(14,540)
## **GROW** @depth 3: [8,0.545455], n=(496,206)
## r=1000 d=[0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0]; n=(83,289,496,206,441,113,14,540,165,366,1518,812,748,472,1068,289)
## **PRUNE** @depth 5: [3,0]
## **GROW** @depth 4: [9,0], n=(146,106)
## **CPRUNE** @depth 1: var=9, val=0.333333->0, n=(521,875)
## **PRUNE** @depth 3: [6,0.46]
## **GROW** @depth 4: [2,0], n=(449,1069)
## **PRUNE** @depth 3: [2,0]
## **GROW** @depth 3: [5,0.155116], n=(126,668)
## **CPRUNE** @depth 0: var=5, val=0.327529->0.344828, n=(1804,5816)
## **GROW** @depth 5: [5,0.551724], n=(206,226)
## **PRUNE** @depth 5: [5,0.551724]
## **GROW** @depth 4: [9,0.666667], n=(290,76)
## **PRUNE** @depth 4: [9,0.666667]
## **GROW** @depth 4: [5,0.293047], n=(476,410)
## **GROW** @depth 5: [6,0.5], n=(218,258)
## r=2000 d=[0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0]; n=(518,323,154,218,258,418,534,446,165,367,366,910,1114,502,1038,289)
##
## Sampling @ nn=0 pred locs:
## **PRUNE** @depth 5: [6,0.5]
## **GROW** @depth 5: [9,0.333333], n=(79,75)
## **GROW** @depth 5: [3,0], n=(560,554)
## **PRUNE** @depth 5: [3,0]
## **GROW** @depth 4: [5,0.793103], n=(117,48)
## **GROW** @depth 5: [6,0.43], n=(149,327)
## **PRUNE** @depth 4: [2,0]
## **GROW** @depth 4: [9,0.333333], n=(62,305)
## **GROW** @depth 4: [2,0], n=(337,832)
## r=1000 d=[0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0]; mh=7 n=(576,378,102,52,147,329,572,716,380,123,51,76,330,337,832,854,456,1038,271)
## **PRUNE** @depth 5: [5,0.172414]
## **GROW** @depth 6: [5,0.655172], n=(383,655)
## **PRUNE** @depth 6: [5,0.655172]
## **GROW** @depth 6: [5,0.172414], n=(76,253)
## **GROW** @depth 3: [2,0], n=(266,450)
## **GROW** @depth 5: [9,0.666667], n=(87,29)
## r=2000 d=[0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0]; mh=8 n=(576,378,199,76,253,171,503,266,450,379,87,29,58,76,331,337,832,998,413,937,271)
## **PRUNE** @depth 5: [9,0.666667]
## **GROW** @depth 7: [2,0], n=(23,53)
## **PRUNE** @depth 3: [2,0]
## **PRUNE** @depth 4: [2,0]
## **GROW** @depth 3: [9,0.333333], n=(97,282)
## **PRUNE** @depth 7: [5,0.137931]
## **GROW** @depth 5: [6,0.79], n=(696,302)
## **GROW** @depth 6: [8,0.545455], n=(339,357)
## **GROW** @depth 4: [6,0.49], n=(214,289)
## **GROW** @depth 4: [2,0], n=(336,828)
## **PRUNE** @depth 4: [6,0.49]
## r=3000 d=[0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0]; mh=8 n=(576,378,199,103,226,171,503,716,30,377,114,71,83,348,336,828,362,377,328,357,880,257)
## **PRUNE** @depth 4: [2,0]
## **GROW** @depth 7: [6,0.83], n=(82,21)
## **GROW** @depth 6: [6,0.69], n=(118,53)
## **GROW** @depth 3: [6,0.48], n=(186,530)
## **PRUNE** @depth 3: [6,0.48]
## **PRUNE** @depth 6: [6,0.7]
## **PRUNE** @depth 7: [6,0.81]
## r=4000 d=[0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0]; mh=8 n=(576,378,199,103,226,171,503,716,30,376,111,74,83,349,1164,606,346,399,272,681,257)
## **GROW** @depth 4: [2,0], n=(337,832)
## **PRUNE** @depth 4: [2,0]
## **GROW** @depth 3: [2,0], n=(266,450)
## **GROW** @depth 6: [5,0.155116], n=(32,139)
## **GROW** @depth 6: [9,0.666667], n=(200,144)
## **PRUNE** @depth 6: [5,0.155116]
## r=5000 d=[0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0]; mh=8 n=(576,378,199,103,226,171,503,266,450,26,327,121,58,77,351,1169,549,481,200,144,293,681,271)
## Grow: 12.95%, Prune: 6.928%, Change: 22.75%, Swap: 8.26%
bike_pred_bcart1 <- predict(bike_mod_bcart1,
task = bike_task,
subset = bike_train_task1)
bike_pred_bcart1
## Prediction: 7620 observations
## predict.type: response
## threshold:
## time: 84.98
## id truth response
## 9165 9165 40 38.4029
## 652 652 122 158.7389
## 8110 8110 37 38.4029
## 3494 3494 153 270.1145
## 2278 2278 379 385.9544
## 8484 8484 356 270.1145
## ... (#rows: 7620, #cols: 3)
Расчет ошибок моделей:
reg_ms1 <- list(rmse, mape, mae)
performance(bike_pred_bcart1,
measures = reg_ms) %>%
round(2)
## rmse mape mae
## 128.69 1.53 89.90
Вывод: как оказалось, из-за технических особенностей наших компютеров, наша команда не может рассчитать модели и их ошибки из-за очень длительного процесса вычиления. Спрогнозировав 3 модели видно, что в целом цедевая ошибка Mape гораздо выше, чем у моделей рассчитаных вручную. К сожалению, пришлось ограничится только этими моделями.
Прогноз на обучающем периоде:
model_2_pred <- bike_train %>%
add_predictions(model_2, var = "Predicted") %>%
dplyr::select(count, Predicted)%>%
mutate(Predicted = Predicted^(1/0.274))
model_3_pred <- bike_train %>%
add_predictions(model_3, var = "Predicted") %>%
dplyr::select(count, Predicted)
model_4_pred <- bike_train %>%
add_predictions(model_4, var = "Predicted") %>%
dplyr::select(count, Predicted)
model_5_pred <- bike_train %>%
add_predictions(model_5, var = "Predicted") %>%
dplyr::select(count, Predicted)
model_6_pred <- bike_train %>%
add_predictions(model_6, var = "Predicted") %>%
dplyr::select(count, Predicted)
model_7_pred <- bike_train %>%
add_predictions(model_6b, var = "Predicted") %>%
dplyr::select(count, Predicted) %>%
mutate(Predicted = Predicted^(1/0.308))
model_8_pred <- bike_train %>%
add_predictions(model_6b1, var = "Predicted") %>%
dplyr::select(count, Predicted) %>%
mutate(Predicted = Predicted^(1/0.308))
model_9_pred <- bike_train %>%
add_predictions(model_6b2m, var = "Predicted") %>%
dplyr::select(count, Predicted) %>%
mutate(Predicted = Predicted^(1/0.314))
model_rtree <- bike_train %>%
add_predictions(rtree, var = "Predicted") %>%
dplyr::select(count, Predicted)
model_rtree_cp <- bike_train %>%
add_predictions(rtree_overfit1, var = "Predicted") %>%
dplyr::select(count, Predicted)
model_rctree <- bike_train %>%
add_predictions(rctree, var = "Predicted") %>%
dplyr::select(count, Predicted)
Расчёт ошибок прогнозов, построенных на основе отобранных ранее ручных и “деревянных” моделей:
errors <- function(actual, predicted) {
error <- actual - predicted
abs_error <- abs(actual - predicted)
percent_error <- abs_error/actual
ME <- mean(error, na.rm = T)
MAE <- mean(abs_error, na.rm = T)
MAPE <- mean(percent_error, na.rm = T)
MSE <- mean(error^2, na.rm = T)
RMSE <- sqrt(MSE)
BIAS <- mean(sum(error)/sum(actual))
tibble(ME, MAE, MAPE, MSE, RMSE, BIAS)
}
models_row <-
bind_rows("2" = model_2_pred, ##ручная
"3" = model_3_pred, ##ручная
"4" = model_4_pred, ##ручная
"5" = model_5_pred, ##ручная
"6" = model_6_pred, ##ручная
"7" = model_7_pred, ##ручная
"8" = model_8_pred, ##ручная
"9" = model_9_pred, ##ручная
"10" = model_rtree, ##дерево решений
"11" = model_rtree_cp, ##дерево решений
"12" = model_rctree, ##дерево решений
.id = "Number")
models_row %>%
group_by(Number) %>%
summarise(errors = list(errors(count, Predicted))) %>%
unnest() %>%
arrange(MAPE) %>%
knitr::kable()
| Number | ME | MAE | MAPE | MSE | RMSE | BIAS |
|---|---|---|---|---|---|---|
| 9 | 33.9135 | 89.0040 | 1.00187 | 19247.07 | 138.734 | 0.177432 |
| 7 | 34.6628 | 89.0989 | 1.00367 | 19335.14 | 139.051 | 0.181352 |
| 8 | 34.6190 | 89.2499 | 1.00425 | 19381.08 | 139.216 | 0.181122 |
| 2 | 26.7644 | 91.4387 | 1.12310 | 19418.70 | 139.351 | 0.140028 |
| 11 | 0.0000 | 63.0498 | 1.30552 | 7971.81 | 89.285 | 0.000000 |
| 12 | 0.0000 | 82.6160 | 1.34796 | 14562.13 | 120.674 | 0.000000 |
| 10 | 0.0000 | 93.9165 | 1.73647 | 17794.19 | 133.395 | 0.000000 |
| 6 | 0.0000 | 107.6542 | 3.00664 | 21329.52 | 146.046 | 0.000000 |
| 5 | 0.0000 | 107.8052 | 3.07060 | 21352.95 | 146.126 | 0.000000 |
| 4 | 0.0000 | 107.7683 | 3.07074 | 21346.71 | 146.105 | 0.000000 |
| 3 | 0.0000 | 107.8035 | 3.07125 | 21341.75 | 146.088 | 0.000000 |
Рассчёт ошибок автоматических моделей
models_row1 <-
bind_rows("randomForest" = performance(bike_pred_lm1,
measures = reg_ms) ,
"lm" = performance(bike_pred_lm,
measures = reg_ms),
"bcart" = performance(bike_pred_bcart1,
measures = reg_ms),
.id = "Number")
models_row1 %>%
group_by(Number) %>%
unnest() %>%
arrange(mape) %>%
knitr::kable()
| Number | rmse | mape | mae |
|---|---|---|---|
| randomForest | 66.5993 | 0.798577 | 44.5661 |
| bcart | 128.6891 | 1.528816 | 89.8989 |
| lm | 144.9483 | 3.160999 | 107.4409 |
Вывод:
Таким образом, можно выделить “модель-победительницу”на обучающей выборке - это автоматическая модель “randomForest”. Она обладает наименьшей MAPE, MSE, RMSE и MAE по сравнению с остальными отобранными моделями.
Второе место берёт ручная модель №9, которая представляет из себя модифицированную модель №6, проверенной на выборке без значимых выбросов и лишённой гетеродоксичносности. Она показала наименьшую ошибку MAPE c небольшим отрывом от остальных моделей и наименьшие знечения остальных ошибок среди остальных оставшихся моделей.
Третье место забирает модель №7 - это модель №6 с переменой atemp вместо temp. Она обладает чуть большей MAPE, MAE, MSE, RMSE, чем у первых двух моделей в рейтинге, но чуть меньшими ME и BIAS, чем у модели, занявшей второе место.
Остальные модели являются “аутсайдерами”.
Построение прогнозов выбранными моделями на тестовой выборке:
model_test_7 <- lm(count^0.308 ~ temp + daytime + I(daytime^2)+ workingday + I(month^2)+ humidity + daytime*workingday + month + weather*daytime, data = bike_test)
model_test7_pred <- bike_test %>%
add_predictions(model_test_7, var = "Predicted") %>%
mutate(Predicted = Predicted^(1/0.308)) %>%
dplyr::select(count, Predicted)
model_test_9 <- lm(count^0.314 ~ atemp + daytime + I(daytime^2) + humidity + daytime*workingday + month*daytime + weather*daytime + humidity*month + humidity*daytime, data = bike_test)
model_test9_pred <- bike_test %>%
add_predictions(model_test_9, var = "Predicted") %>%
mutate(Predicted = Predicted^(1/0.314)) %>%
dplyr::select(count, Predicted)
Расчёт ошибок прогнозов моделей №7 и №8 на тестовой выборке:
errors <- function(actual, predicted) {
error <- actual - predicted
abs_error <- abs(actual - predicted)
percent_error <- abs_error/actual
ME <- mean(error, na.rm = T)
MAE <- mean(abs_error, na.rm = T)
MAPE <- mean(percent_error, na.rm = T)
MSE <- mean(error^2, na.rm = T)
RMSE <- sqrt(MSE)
BIAS <- mean(sum(error)/sum(actual))
tibble(ME, MAE, MAPE, MSE, RMSE, BIAS)
}
models_row_test <-
bind_rows(
"7" = model_test7_pred,
"9" = model_test9_pred,
.id = "Number")
models_row_test %>%
group_by(Number) %>%
summarise(errors = list(errors(count, Predicted))) %>%
unnest() %>%
arrange(MAPE) %>%
knitr::kable()
| Number | ME | MAE | MAPE | MSE | RMSE | BIAS |
|---|---|---|---|---|---|---|
| 7 | 23.6644 | 89.3452 | 1.05954 | 18234.0 | 135.033 | 0.122407 |
| 9 | 22.8688 | 89.4207 | 1.07390 | 18355.1 | 135.481 | 0.118291 |
Прогноз “randomForest” на тестовой выборке:
bike_mod_lm11 <- train(lrn_lm1, task = bike_task,
subset = bike_train_task1)
bike_pred_lm11 <- predict(bike_mod_lm11,
task = bike_task,
subset = bike_test_task1)
reg_ms <- list(rmse, mape, mae)
performance(bike_pred_lm11,
measures = reg_ms) %>%
round(2)
## rmse mape mae
## 111.53 1.31 74.86
Вывод: наилучшей регрессионной моделью на тестовой выброке стала модель №7. Она обладет наименьшей MAPE, совсем незначительно большими значениями других ошибок по сравнению с моделью, занявшей второе место.второе место - модель №9, а третье - автоматическая модель “randomForest”.
Таким образом, удалось определить наилучшую регрессионную модель для прогнозирования, которая показала наилучшие показатели ошибок на тестовой выброке - это модель №9, которая включала в себя множественную регрессию переменных, которая была протестирована на выборке, лишённой значимых выбросов, взаимосвязь переменных, а также избавленной от возможной гетеродокичности. Эта модель заняла третье место по точности прогнозирования на обучающей выборке, но это не помешало ей занять первое место на тестовой выборке. В заключение хотелось бы скзать, что найденная нами регрессионная модель позволяет довольно точно предсказывать значения целевой переменной - число актов аренды велосипедов в зависимости от некоторых внешних факторов, которые определены другими переменными в датасете.